Данные об университетах.
Признаки:
Выберем некоторые из этих признаков.
uni <- read.csv2("/Users/Vladimir/Dropbox/Statistics_Autmn_2016/PCA/I.csv")
uni <- uni[,c("X", "PPIND", "NEW10", "R_B_COST","ADD_FEE", "BOOK", "SAL_ALL", "NUM_ALL" ,"APP_ACC", "SF_RATIO", "IN_STATE", "OUT_STAT")]
uni <- na.omit(uni)
nums <- c("NEW10", "R_B_COST","ADD_FEE", "BOOK", "SAL_ALL", "NUM_ALL" ,"APP_ACC", "SF_RATIO", "IN_STATE", "OUT_STAT")
uni[,nums] <- scale(log(uni[,nums]))
Матриксплот:
pairs.panels(uni[,nums],
pch=21,
lm=TRUE,
lwd = 1,
ellipses = FALSE)
По скаттерплотам можно предположить, что существует два кластера. Это хорошо видно, например, на графике NUM_ALL vs. IN_STATE.
Модель – смесь нормальных распределений: \[ p(x) = \sum_{k = 1}^G w_k p_k(x), ~p_k(x) = \phi(x; \mu_k, \Sigma_k). \]
Хотим найти параметры максимизируя функцию правдоподобия: \[ \mathcal{L}(\theta; \mathbf{x}) = \prod_{i = 1}^n \sum_{k = 1}^G w_k \phi(x_i; \mu_k, \Sigma_k). \]
Подробно разберем случай смеси двух распределений.
Смесь: \[ p(x) = (1 - \pi)\phi(x; \mu_1, \Sigma_1) + \pi\phi(x; \mu_2, \Sigma_2). \]
Логарифм функции правдоподобия:
\[ \mathcal{l}(\theta; \mathbf{x}) = \sum_{i = 1}^n log((1 - \pi)\phi(x_i; \mu_1, \Sigma_1) + \pi\phi(x_i; \mu_2, \Sigma_2)). \]
Численно максимизировать \(\mathcal{l}(\theta; x)\) сложно, поэтому вводятся ненаблюдаемые переменные \(\Delta_i\), принимающие значения \(0\) или \(1\). Если \(\Delta_i = 1\), то элемент выборки \(x_i\) пришел к нам из первого распределения, иначе – из второго. Логарифм функции правдопадобия перепишется так:
\[ \mathcal{l}(\theta; \mathbf{x}, \mathbf{\Delta}) = \sum_{i = 1}^n (1 - \Delta_i) \log \phi(x_i; \mu_1, \Sigma_1) + \Delta_i \log \phi(x_i; \mu_2, \Sigma_2) + \sum_{i = 1}^n (1 - \Delta_i) \log(1 - \pi) + \Delta_i \log \pi. \]
Так как значения \(\Delta_i\) неизвестны, заменим их на условные мат. ожидания:
\[ \gamma_i(\theta) = \mathbb{E}(\Delta_i|\theta, \mathbf{x}). \]
Далее воспользуемся ЕМ-алгоритмом, чтобы найти оценки параметров.
Пусть у нас есть начальные оценки параметров: \(\hat{\mu}_1\), \(\widehat{\Sigma}_1\), \(\hat{\mu}_2\), \(\widehat{\Sigma}_2\), \(\hat{\pi}\).
Expectation Step
\[ \hat{\gamma}_i = \frac{\hat{\pi}\phi(x_i; \mu_2, \Sigma_2)}{(1 - \pi)\phi(x_i; \mu_1, \Sigma_1) + \pi\phi(x_i; \mu_2, \Sigma_2)}. \]
Maximization Step
\[ \hat{\mu}_1 = \frac{\sum_{i = 1}^n (1 - \hat{\gamma}_i)x_i}{\sum_{i = 1}^n (1 - \hat{\gamma}_i)}, \] \[ \hat{\mu}_2 = \frac{\sum_{i = 1}^n \hat{\gamma}_i x_i}{\sum_{i = 1}^n \hat{\gamma}_i}. \] Если не делать никаких предположений о виде ковариационных матриц, то их оценки: \[ \widehat{\Sigma}_1 = \frac{\sum_{i = 1}^n (1 - \hat{\gamma}_i) (x_i - \hat{\mu}_1)(x_i - \hat{\mu}_1)^{\mathrm{T}}}{\sum_{i = 1}^n (1 - \hat{\gamma}_i)}, \] \[ \widehat{\Sigma}_2 = \frac{\sum_{i = 1}^n \hat{\gamma}_i (x_i - \hat{\mu}_1)(x_i - \hat{\mu}_1)^{\mathrm{T}}}{\sum_{i = 1}^n \hat{\gamma}_i}, \] \[ \hat{\pi} = \frac{1}{n}\sum_{i = 1}^n \hat{\gamma}_i. \]
Продолжаем пока не сойдемся.
Ковариационные матрицы определяют вид эллипсов.
Найдем спектральное разложение ковариационной матрицы: \[ \Sigma_k = \lambda_k D_k A_k D_k^\mathrm{T}, \] где \(D_k\) – матрица, составленная из собственных векторов \(\Sigma_k\), \(A_k\) – диагональная матрица, компоненты которой пропорциональны собственным числам, \(\lambda_k\) – коэффициент пропорциональности.
Такое разложение можно интерпретировать геометрически: \(D_k\) отвечает за ориентацию эллипса, \(A_k\) – за его форму, \(\lambda_k\) – за его объем. Фиксируя различные признаки эллипсов, получаем различные модели. В случае размерности большей 2 модель обозначают тремя буквами, например, EVI означает, что у эллипсов одинаковый объем, возможно разные формы, а корреляции равны 0.
Байесовский информационный критерий.
Идея: давайте штрафовать за большое число параметров.
\[ BIC = -2 \hat{\mathcal{l}} + k \log(n), \] где \(n\) – число точек, \(k\) – число параметров
raw.data <- uni[,nums]
Посмотрим на байесовский информационный критерий.
BIC <- mclustBIC(raw.data)
plot(BIC)
summary(BIC)
## Best BIC values:
## EEE,2 EEE,5 EEE,3
## BIC -2899.265 -2908.343604 -2913.00935
## BIC diff 0.000 -9.078216 -13.74397
BIC рекомендует предположить, что у нас два кластера, а ковариационные матрицы одинаковы. Но по матриксплотам выше видно, что точнее будет взять VEV, так как объемы разные, формы одинаковые, корреляции тоже разные (например, NEW_10 vs. R_B_COST).
m <- Mclust(raw.data, G = 2, modelNames = "VEV")
Результат классификации.
plot(m, "classification")
clusplot(raw.data, m$classification, main='2D representation of the Cluster solution',
color=TRUE, shade=TRUE,
labels=2, lines=0)
Идея – для каждой точки определяем насколько она похожа на остальные точки из своего кастера по сравнению с точками из других кластеров. Например, в смысле евклидовой метрики.
Пусть \(a(i)\) – среднее расстояние от \(i\)-й точки до остальных в кластере (cohesion), \(b(i)\) – самое маленькое из средних расстояний до точек из других кластеров (separation).
Силуэт определяется так:
\[ s(i) = \frac{b(i) - a(i)}{max\{a(i), b(i)\}}. \]
Чем ближе \(s(i)\) к \(1\), тем более уверенно можно сказать, что точка в нужном кластере. Среднее значение \(s\) – мера того, как хорошо кластеризовали.
dissE <- daisy(raw.data)
dE2 <- dissE^2
sk2 <- silhouette(m$classification, dE2)
plot(sk2)
Кроме того, силуэт можно использовать как один из способов предположить количество кластеров. Если выбрали слишком большое или слишком маленькое k, некоторые селуэты окажутся заметно уже чем остальные. Например, если возьмем 6 кластеров для наших данных:
m6 <- Mclust(raw.data, G = 6, modelNames = "EEE")
dissE <- daisy(raw.data) #считаем среднее расстояние от каждой точки до остальных (внутри и вне кластера)
dE2 <- dissE^2
sk2 <- silhouette(m6$classification, dE2)
plot(sk2)
cl <- kmeans(raw.data, 2)
plot(raw.data, col = c("red", "blue")[cl$cluster])
Картинка в главных компонентах:
cl <- kmeans(scale(raw.data), 2)
clusplot(scale(raw.data), cl$cluster, main='2D representation of the Cluster solution',
color=TRUE, shade=TRUE,
labels=2, lines=0)
Silhouette plot
dissE <- daisy(raw.data)
dE2 <- dissE^2
sk2 <- silhouette(cl$cluster, dE2)
plot(sk2)
Если верить силуэту, данные разделились чуть хуже чем при model-based clustering (среднее s(i) 0.5 против 0.51 в случае model-based).
На каждом шаге считается расстояние до “дальнего соседа”.
\[ R(W, S) = \max_{s \in S, w \in W} \rho(w,s) \]
#euclid dist
d <- dist(raw.data)
hc.cl <- hclust(d, method = "complete")
plot(hc.cl)
plot(raw.data, col = c("red", "blue")[cutree(hc.cl, 2)])
Здесь расстояние – то, насколько увеличится сумма квадратов расстояний до центра при слиянии двух кластеров, цена слияния.
\[ \Delta(A, B) = \sum_{i \in A \cup B}\|x_i - \mu_{A \cup B} \|^2 - \sum_{i \in A}\|x_i - \mu_{A} \|^2 - \sum_{i \in B}\|x_i - \mu_{B} \|^2 \]
#euclid dist
d <- dist(raw.data)
hc.ward <- hclust(d, method = "ward.D2")
plot(hc.ward)
Так как метод Варда на каждом шаге выбирает минимальную цену слияния, большие скачки могут означать, что мы объединяем точки не внутри кластеров, а сами кластеры. Тогда по картинке выше можно решить, что у нас два кластера.
plot(raw.data, col = c("red", "blue")[cutree(hc.ward, 2)])
На самом деле университеты делятся на государственные и частные. Выглядит это так:
pairs.panels(uni[,nums],
pch=21 +as.numeric(uni$PPIND),
bg = c("red", "blue")[uni$PPIND],
lm=TRUE,
lwd = 1,
ellipses = FALSE)
Проверим, насколько ошиблись методы выше.
table(uni$PPIND, m$classification)
##
## 1 2
## 1 75 0
## 2 3 42
Всего три ошибки.
table(uni$PPIND, cl$cluster)
##
## 1 2
## 1 2 73
## 2 42 3
Ошибок чуть больше чем в model-based. Если вернуться к матриксплоту выше, видно, что model-based справился, а k-means – не совсем.
table(uni$PPIND, cutree(hc.cl, 2))
##
## 1 2
## 1 67 8
## 2 2 43
table(uni$PPIND, cutree(hc.cl, 2))
##
## 1 2
## 1 67 8
## 2 2 43
Количество ошибок больше чем в двух методах выше. По матриксплотам выше видно, что в случае иерархической кластеризации получилось что-то похожее на k-means.